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Abstract 

Biochemical reactions in living systems occur in complex, heterogeneous media with total concentrations of 
macromolecules in the range of 50 - 400 Molecular species occupy a significant fraction of the immersing 
medium, up to 40% of volume. Such compfex and volume-occupied environments are generally termed 'crowded' 
and/or 'confined'. In crowded conditions non-specific interactions between macromolecules may hinder diffusion - 
a major process determining metabolism, transport, and signaling. Also, the crowded media can alter, both 
qualitatively and quantitatively, the reactions in vivo in comparison with their in vitro counterparts. This review 
focuses on recent developments in particle-based Brownian dynamics algorithms, their applications to model 
diffusive transport in crowded systems, and their abilities to reproduce and predict the behavior of 
macromolecules under in vivo conditions. 



Introduction 

Intracellular organelles are packed with small solutes, 
macromolecules, membranes and skeletal proteins. Such 
crowding is characteristic not only of the cell's interior but 
also for extracellular tissues [1]. The kinetics and thermo- 
dynamics of macromolecular processes and biochemical 
reactions taking place in vivo are known to be affected by 
complex, volume-occupied environments [2-5]. 

More and more experimental methods have become 
available to investigate macromolecules under in vivo 
conditions. These include nuclear magnetic resonance 
[6-8] and fluorescence spectroscopies [9]. Techniques 
such as SPT (single-particle tracking) [10,11], FRAP 
(fluorescence recovery after photobleaching) [12,13], FCS 
(fluorescence correlation spectroscopy) [14] have been 
applied to measure diffusion constants of macromole- 
cules in the cytoplasm and membranes. All experiments 
show that diffusion of proteins in vivo is significantly 
reduced compared to dilute conditions. In the cytoplasm 
of eukaryotic cells, diffusion of both large and small 
molecules is slowed down three to four times [1]. FRAP 
measurements of the diffusion coefficient of GFP in the 
Escherichia coli (E. coli) cytoplasm [15,16] yield about 
10 times smaller values than at infinite dilution in water. 
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Theories of diffusion for colloidal soft matter are well 
established. Within their framework, the behavior of col- 
loidal systems (both dilute and concentrated), consisting 
of mesoscopically large colloidal particles dispersed in 
low-molecular-weight solvent, can be simulated, repro- 
duced, and more importantly predicted [17-21]. In cell 
biology, however, an appropriate theory of diffusion phe- 
nomena is still lacking. Complications arise due to the 
cellular heterogeneity. Moreover, there are different cell 
types whose compositions vary over space and time on 
the scales associated with diffusive motion and the 
knowledge of the surroundings of diffusing species is 
essential for successful applications of theoretical 
approaches. Unfortunately, even for the best-studied sys- 
tems, experimental data is too scarce to describe elemen- 
tary processes involving macromolecules, their timing 
and spatial distribution. 

Simulation methods used to study the in vivo diffusion 
of macromolecules should take into account their bio- 
logical localization, sizes and shapes, and positions in 
three dimensional space. Intermolecular interactions 
(both specific and non-specific, repulsive and attractive) 
should be also included [22] . Ideally, the modeling should 
be performed at atomic resolution. However, due to the 
[im size of an average cell and a large number of compo- 
nents occupying its volume, reductionist (coarse-grained) 
approaches are unavoidable. The simulation methods 
[22] that can be applied to directly account for the 
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crowded conditions, sizes and shapes of molecules and 
their interactions, as well as other features of the cellular 
environment, include different variants of particle-based 
Brownian dynamics (BD) [23-26]. In this review, we 
describe recent applications of BD to study diffusion in 
the crowded, in-vivo-like systems. Examples of different 
modeling approaches, either coarse-grained or more ela- 
borate, are presented and discussed. 

Diffusion in dilute and crowded solutions 

In dilute solutions, diffusional properties of molecules 
treated as rigid bodies are determined by their sizes, 
shapes, temperature, and solvent viscosity. The complete 
information required to describe translational and rota- 
tional motions of Brownian particles in dilute solutions 
is contained in their diffusion tensors. Single-particle 
diffusion tensor, D, is represented with a symmetric, 
6x6 matrix containing 3x3 blocks related to the transla- 
tional (tt) and rotational diffusivities (rr) and their cou- 
plings (rt, tr) [26]: 



v J 

The diffusion tensors of molecules can be obtained 
from rigid-body hydrodynamic calculations [27] . Transla- 
tional diffusion coefficient is an average over the diagonal 
terms of D tt . The rotational dynamics of a free diffusing 
molecule can be predicted from the eigenvalues of D rr 
[28]. For a spherical molecule, diffusion is isotropic, the 
diffusion tensor is diagonal, and the molecule can be 
described with a single translational and a single rota- 
tional diffusion coefficient. Translational diffusion coeffi- 
cient (£) t,aws ) determines the time dependence of the 
diffusion of the particle's center of mass. The mean 
square displacement (MSD) of a freely diffusing (in three 
dimensions) particle's center of mass changes linearly 
with time t according to the following equation [29]: 

MSD(t) = 6 ■ D trans ■ t (2) 

with the D am as a proportionality constant. 

In biological setup, solutes interact with their environ- 
ments and other molecules which affects their Brownian 
motion, thus single-particle diffusion tensors are no 
longer the unique determinants of diffusion. Neverthe- 
less, the relation between the MSD and time t (not 
necessarily linear) can still be used to measure the 
apparent translational diffusion of molecules [4,11]. 

Particle-based Brownian dynamics algorithms 

Brownian dynamics is a stochastic simulation approach 
with continuum space and time. In BD the molecules 



exhibit noise as they are propagated according to the 
Langevin equation [30]. The equation contains random 
and frictional forces that are intended to represent the 
interaction between the diffusing molecules and solvent 
(treated implicitly). High-friction limit is assumed; 
momentum relaxation is much faster than position relaxa- 
tion, thus the equations of motions do not contain parti- 
cles' velocities. Below, we briefly describe two common 
BD algorithms. 

Brownian dynamics of rigid, arbitrarily shaped objects 

The propagation scheme for a number of rigid bodies, 
described either at atomic level or with coarse-grained 
models, can be written as [23,26]: 

x i =x° i +-^-D i M l +R l (M) (3) 

where i runs over molecules, At is the time step, and 
x is the vector describing the position of the center of 
diffusion ( r) and orientation (^) of the i-th molecule, 
M- M is a generalized force acting on a given particle 
having two components: the force acting on the diffu- 
sion center ( p ) and the torque (f\,M = (F,f \ . ^ (At) is 
a random displacement vector resulting from the Brow- 
nian noise, with zero mean and the variance-covariance 
given with: 

< Rj ( At ) Rj ( At ) >= 2D i S i) A£ (4) 

Here, D is the precomputed and constant in the simu- 
lation 6x6 diffusion tensor of the 2-th molecule at its 
diffusion center. Random displacements of particles can 
be computed via Cholesky decomposition of D [23]. 
Intermolecular hydrodynamic interactions between dif- 
fusing bodies are neglected. 

A simplified version of the above algorithm [25], 
which assumes isotropic properties of molecules and 
replaces their diffusion tensors with translational and 
rotational diffusion coefficients, has been widely used. 
Apart from treating the diffusing molecules as spheres, 
this simplification also ignores the coupling between 
translational and rotational diffusion. It is however ques- 
tionable, whether using pre-averaged diffusion tensors is 
justified for nonglobular molecules (e. g., see the discus- 
sion on the influence of anisotropy on transport and 
association rates of model molecules [31] or on the dif- 
fusion of an ellipsoid near planar surfaces [27]). 

Brownian dynamics with hydrodynamic interactions 

When hydrodynamic interactions (HI) are considered in 
BD simulations, the diffusion tensor of the entire system 
(6Nx6N, with N being the number of spherical mole- 
cules or the number of spherical pseudo-atoms in case 
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of more elaborate coarse-grained molecular models [32]) 
is used to propagate particles. This tensor depends on 
the system's configuration and varies with time. The 
propagation scheme can be written as [23,24]: 



x = x° +(VD)At + — DM + R(At) 



(5) 



When the Rotne-Prager-Yamakwa [33,34] form of the 
diffusion tensor is used (see below) the gradient term in 
the above equation (VD) vanishes. The 6N-dimensional 
random displacement vector has zero mean and fulfills 
< (At) (At) > = 2D At. Again, random displacements 

of particles can be computed via the Cholesky decompo- 
sition of D [23]. However, the decomposition of D must 
be performed at each BD step due to its dependence on 
the positions of particles. 

Intermolecular interactions 

The influence of environments on the diffusion of 
macromolecules manifests through intermolecular forces 
acting on individual particles. The types and functional 
forms of interactions included in BD simulations vary 
depending on the level of detail used to describe the 
studied systems. Different types of intermolecular inter- 
actions that can be modeled in BD simulations are 
briefly described below. 

Electrostatics 

When diffusing molecules are described in a reduced 
way with spherical subunits (e. g., as in case of model 
spherical proteins [35-37] or bead-models of flexible 
polymers [38]), the DLVO (Derjaguin - Landau - 
Verwey - Overbeek) approach (or the Debye-Hiickel 
approximation), commonly applied in the studies of col- 
loid suspensions, can be used [39]. Spherical molecules 
(of radii R it Rj) with central net charges (Q it Qj) interact 
via pairwise additive potentials: 



V 



QiQi 



Ane 0 e ( 1 + feRj ) ( 1 + kRj \ r { 



-K<r,:-R,-R:) 

e >' (6) 



Where e a is the vacuum permittivity, e is the dielectric 
constant of the immersing medium, k is the inverse of 
the Debye screening length, and is the separation 
between molecules. However, the above description is 
not valid for highly charged systems such as nucleic 
acids where the linearized form of the Poisson-Boltz- 
mann equation [40] is not applicable and there is strong 
electrostatic coupling between the molecular and ionic 
species. In such cases, the net charges can be replaced 
with renormalized charges [41]. 

A more sophisticated approach to treat electrostatic 
interactions in BD simulations of rigid molecules with 



atomic detail [42,43] is based on the model of effective 
charges developed by Gabdoulline and Wade [44]. In 
this approach, each molecule is described with a limited 
set of optimized, discrete, effective charges immersed in 
an uniform dielectric. These effective charges are 
derived by fitting the electrostatic potential resulting 
from the Debye-Hiickel approximation to the external 
molecular potential obtained from the numerical solu- 
tion of the Poisson-Boltzmann equation. Electrostatic 
energies (and thus forces) are evaluated using the values 
of the electrostatic potential (derived from the solution 
of the Poisson-Boltzmann equation on a three-dimen- 
sional grid) generated by one molecule at positions of 
effective charges of the second molecule. The electro- 
static potential grids and effective charges are precom- 
puted so when multiple molecules are simulated, the 
electrostatic potential grids need to be translated and 
reorientated at each simulation step along with translat- 
ing and rotating molecules. This approach ignores the 
dielectric polarization effects that occur when two mole- 
cules come close to each other. These effects may be 
particularly important in the crowded systems, where 
intermolecular distances are comparable with the sizes 
of molecules. Polarization can be accounted for by using 
an analytical formula proposed by Gabdoulline and 
Wade [44]. The performance of the effective charges 
approach is worsened for highly charged molecules 
in cases where the nonlinear form of the Poisson- 
Boltzmann equation should be used to evaluate electro- 
static potentials. The Debye-Hiickel approximation is 
not valid in these cases, especially in the vicinity of the 
electrostatic field sources [45]. 

Lennard-Jones potentials 

Non-specific interactions between molecules can be 
modeled with the Lennard-Jones type functions that are 
commonly used in molecular dynamics simulations to 
compute the van der Waals interactions between atoms; 
at large intermolecular separations molecules attract 
each other, while at small separations the interactions 
are repulsive and molecules are impenetrable. The 
Lennard-Jones functional form and parameters for BD 
simulations with coarse-grained spherical models [35,37] 
can be, for instance, taken from the ones derived for 
large colloidal spheres [46,47]. 

In BD simulations employing all-atom models, the 
forces acting between the surface atoms of molecules 
can be evaluated using standard (6/12) Lennard-Jones 
potentials: 



vl2 



(7) 
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However, the well depth e L j needs to be treated as a 
free parameter and adjusted to reproduce experimental 
data [42]. 

Hydrophobic effects 

Nonpolar (hydrophobic) interactions can be included in 
atomistic BD simulations by assuming their proportion- 
ality to the amount of the solvent accessible surface area 
(SASA) that gets buried when molecules come close to 
each other. However, computing of SASA can be slow 
and to make SASA-based approaches efficient, approxi- 
mations are needed. Examples of how to incorporate the 
hydrophobic effects based on SASA in BD simulations 
can be found in the work of Elcock and McCammon 
[48] or Gabdoulline and Wade [49]. 

While SASA-based models are commonly used to 
evaluate nonpolar interactions, it has been shown that 
including the solvent accessible volume (SAV) terms 
may be essential to accurately model the nonpolar solva- 
tion, especially at atomic-length scales [50]. Accurate 
treatment of nonpolar energies and forces is crucial for 
understanding any biomolecular process that involves 
changes in solvent accessibility. The described inability 
of SASA-only models [50] to reproduce and predict sol- 
vation energies and forces indicates that there is a need 
to include a more complete theory of nonpolar solvation 
in atomistic BD simulations. 

It was also proposed that hydrophobic interactions can 
be accounted for by modifying the van der Waals inter- 
actions and assigning more favorable interaction ener- 
gies to nonpolar surface atoms than to polar ones [42]. 
With such an approach hydrophobic interactions can be 
evaluated more rapidly, however its correctness is ques- 
tioned (see the discussion in [42]). 

Hydrodynamic interactions 

In BD simulations water molecules are not treated expli- 
citly but the influence of solvent on the diffusion of sus- 
pended molecules can be included by proper treatment of 
hydrodynamic interactions. However, even though their 
importance is widely recognized [32,37,51], accounting for 
HI is no simple matter due to their long-range and many- 
body character [52] resulting in high computational cost. 
Therefore, one is often tempted to neglect HI and, conse- 
quently, the correlations of motions of molecules. 

In the Ermak-McCammon algortihm [23] (Equation 5), 
hydrodynamic interactions between spherical particles 
can be included by using the Oseen [53], Rotne-Prager 
[33] or Rotne-Prager-Yamakawa [34] forms of diffusion 
tensors (without HI the diffusion tensor is represented by 
a diagonal matrix). At each step of a BD simulation, the 
diffusion tensor of the whole system is factorized using 
the Cholesky decomposition [23] to obtain hydrodynami- 
cally correlated displacements of particles. The Cholesky 



decomposition scales as N 3 (with N being the number of 
particles), while evaluating direct intermolecular interac- 
tions, that are assumed pairwise, scales as N 2 . It is possi- 
ble to lower the cost of evaluating HI using the 
Chebyshev approximation proposed by Fixman [54], but 
the Hi-related computational overhead is still consider- 
able (N ). Recently, Geyer and Winter proposed a faster 
algorithm that scales as N and is based on a truncated 
expansion of the hydrodynamic multi-particle correla- 
tions as two-body contributions [55]. 

An additional computational cost arises when HI are 
evaluated in finite simulation boxes containing mole- 
cules. This cost is due to the fact that for evaluating dif- 
fusion tensors Ewald summation has to be used [56,57]. 
While the direct intermolecular interactions in such per- 
iodic systems can be computed using less expensive cut- 
off -based schemes, applying these schemes to evaluate 
HI leads to diffusion tensors that are not positive defi- 
nite (a condition that has to be met to use diffusion 
tensors as covariance matrices [57].) 

The treatment of hydrodynamics within the frame- 
work of the Oseen, Rotne-Prager or Rotne-Prager- 
Yamakwa tensors is appropriate for very dilute systems 
capturing only the far-field, two-body hydrodynamic 
effects. However, in crowded biological systems, the 
many-body HI and lubrication forces (resulting from the 
thin layer of solvent that separates the nearly touching 
particles) may significantly influence diffusion. One 
crude way to include the many-body HI in BD simula- 
tions is to use mean-field approaches [58-60] that are 
based on scaling the diffusion coefficients of particles 
according to local volume fractions, as for example in 
the work of Sun and Weinstein [61]. 

The most advanced approach to treating HI, used 
recently in BD simulations of a biological system [37], 
was developed by Brady and Bossis [17]. Their Stokesian 
dynamics (SD) method includes both far-field and near- 
field many-body contributions to Brownian forces acting 
on particles. Moreover, particular implementations of SD 
enable one to efficiently simulate long-time dynamics of 
large multicomponent systems. For example, the acceler- 
ated SD method, described in [62] scales as N 125 log(N). 

Unfortunately, modeling of HI in BD simulations is 
currently limited to systems composed of coarse-grained 
particles (either spherical molecules or composed of a 
small number of spherical units - pseudo-atoms). While 
such models perform well for colloids, a question yet 
unanswered is whether they are sufficient for crowded, 
heterogeneous biological systems. 

Applications to diffusion in crowded environments 

Typically, the setup of BD simulations aimed at investigat- 
ing diffusion in crowded solutions is the following. 
A number of molecules (usually 10 - 10 ), described 



Diugosz and Trylska BMC Biophysics 201 1, 4:3 
http://www.biomedcentral.eom/2046-1682/4/3 



Page 5 of 9 



either at the atomistic level or with coarse-grained models, 
is contained inside a primary simulation box. Periodic 
boundaries are applied to reduce the edge effects and the 
influence of the finite size of the system. BD trajectories of 
all molecules are generated with the algorithm of choice. 
The BD simulation time scale should be long enough so 
the relative displacements in the system be at least com- 
parable with the dimensions of molecules; the simulation 
time scales range from microseconds to milliseconds, 
while typical time steps are about a few picoseconds. 
When one considers the number of simulated molecules, 
the need to describe them in as detailed way as possible, 
and the complexity of interactions calculated in a single 
simulation step, it becomes clear that BD simulations of 
crowded systems require a lot of computational resources. 

Below we present a few recent (except for one) exam- 
ples of BD simulations of macromolecular diffusion in 
crowded media. We describe the studies that best illus- 
trate the potential and shortcomings of the BD techni- 
que applied to in vivo-like systems and the ones that 
introduce a significant element of novelty. Some of 
these computational studies relate to experiments - 
experimental data were either used for parameterization 
or for direct comparison with simulations. 

Beginnings 

One of the first studies that used the BD technique to 
investigate macromolecular diffusion in congested media 
was that of Dwyer and Bloomfield [63] in 1993. They 
performed BD simulations of probe and self- diffusion of 
concentrated solutions containing short DNA polymers 
(modeled as wormlike chains, i.e. strings of spherical 
beads) and the bovine serum albumin (BSA) protein 
modeled as a sphere. Even though their DNA model 
lacked the detailed all-atom description, it reproduced 
the overall translational and rotational behavior of linear 
DNA molecules [64]. The authors used an algorithm 
similar to that of Ermak-McCammon [23], however, 
they neglected HI, because the pairwise approaches are 
inadequate for concentrated solutions. They also used 
the Debye-Huckel approach to model long-range elec- 
trostatics and the repulsive part of the Lennard-Jones 
potential to ensure volume exclusion. This model, while 
very simple, turned out to be realistic enough to simu- 
late diffusion in the DNA/BSA system and accurately 
reproduce experimental diffusion coefficients at 0.1 M 
ionic strength for the self-diffusion of BSA and the 
probe diffusion of BSA in DNA (both over a wide range 
of BSA/DNA concentrations). 

Simulations of the diffusion of BSA in the DNA solu- 
tion conducted at low ionic strength of 0.01 M gave 
worse agreement with experiments which the authors 
attributed to the lack of HI and limitations of the 
applied electrostatic model. The work of Dwyer and 



Bloomfield is significant because it shows that even very 
simple modeling approaches can accurately predict the 
behavior of rather complex systems. 

We note that the models similar to the one of Dwyer 
and Bloomfield are still used, for example, in computa- 
tional studies of facilitated diffusion of proteins on 
DNA. The description of facilitated diffusion phenom- 
enon and models of non-specific protein-DNA interac- 
tions can be found in references [65,66]. 

Protein-protein association in crowded environments 

Most biological reactions require some sort of diffusion- 
mediated encounter. Here we describe BD applications 
that directly investigate the influence of crowded envir- 
onments on diffusional encounter and association of 
molecules. 

Sun and Weinstein [61] simulated systems of spherical 
particles, mimicking globular proteins. Particles were 
modeled as hard spheres; the authors used the elastic col- 
lision method to effectively deal with particle overlaps 
[67]. Electrostatic interactions were introduced using the 
Debye-Hiickel approach. A notable fact is that the 
authors attempted to model also many-body HI. They 
used the translational part of the Ermak-McCammon 
algorithm with the diagonal form of the diffusion tensor. 
However, they applied an average hydrodynamic correc- 
tion and scaled the diagonal terms of the diffusion tensor 
(i.e., the single particle diffusion coefficients) with func- 
tions that depend on the instantaneous local volume frac- 
tion of each particle [58-61]. The authors investigated 
how the size of background molecules, volume fraction, 
electrostatic interactions, and HI influence the diffusion 
of a tracer molecule. Their study also quantitatively 
described the dependence of diffusion-limited reaction 
rates on crowding in model systems. Using a very 
approximate approach, they showed a significant effect of 
HI on diffusion and association rates. While this review 
focuses on simulations of multimolecular systems, the 
importance of intermolecular HI in modeling the associa- 
tion kinetics in a two-molecular case was also described 
by Frembgen-Kesner and Elcock [32]. A number of sim- 
plifications was used in the described work of Sun and 
Weinstein [61]. Simulated systems were rather homoge- 
neous in comparison with in vivo environments. Mole- 
cules were treated as spheres thus their association was 
non-specific (the model of particles with uniformly reac- 
tive surfaces was used). Anisotropy of particles, which 
may play a role upon association, and specific intermole- 
cular interactions were neglected. However, their work is 
important as it showed the necessity to include HI in BD 
simulations of crowded systems. 

One of the conclusions of the work of Sun and 
Weinstein [61] was that crowding slows down the non- 
specific association of molecules. The two more recent 
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BD studies focused on the influence of crowders on spe- 
cific interactions in protein-protein association and con- 
sidered the fact that partner proteins need to be properly 
oriented to associate. Wieczorek and Zielenkiewicz [68] 
used BD to simulate the diffusional encounter of HEL 
and HyHEL-5 proteins in solutions of hard-sphere crow- 
ders. Both proteins were modeled with atomistic details, 
however, only the excluded volume interactions were 
evaluated; simulation steps leading to overlaps were 
rejected and repeated with different random displace- 
ments, so no systematic forces were explicitly applied. 
The crowders and proteins were modeled as rigid bodies 
with isotropic diffusion tensors (both HEL and HyHEL-5 
were represented as spheres with equivalent hydrody- 
namic radii) and propagated with the algorithm described 
in [25]. Hydrodynamic interactions were neglected. The 
authors concluded that if the definition of an encounter 
complex accounts for the orientation of associating mole- 
cules, crowding can accelerate the association rate. 

Similar observations were made by Kim and Yethiraj 
[69]. Using spherical protein models with non-uniformly 
reactive surfaces, they showed that association of model 
molecules can be either slowed down or accelerated, 
depending on the applied reaction criteria. These 
authors recently used BD to investigate the effects of 
crowding on the association at model membranes [70]. 
Both studies of Kim and Yethiraj [69,70] focus on how 
the excluded volume effects influence the reaction rates. 

In general, the above simulation results are in agree- 
ment with purely theoretical predictions [3]. However, 
they lack a close connection to experiments and the 
ability to reproduce heterogeneous in vivo conditions. 

Modeling bacterial cytoplasm 

In the recent work of Elcock's [43] and Skolnick's [37] 
groups BD was applied to investigate the diffusion in 
cytosol-like systems with the Escherichia coli cytoplasm 
as a model system. The cytoplasm model of McGuffee 
and Elcock [43] contains 50 different types of the most 
abundant proteins and nucleic acids present in the E. coli 
cytoplasm. Additionally, a few molecules of the Green 
Fluorescent Protein (GFP) were included, for the sake of 
comparison with experimental data. All molecules were 
represented with atomistic details. The algorithm of rigid 
bodies was applied in BD but all molecules were treated 
as isotropic, so the average rotational and translational 
diffusion coefficients were used instead of 6x6 diffusion 
tensors [25]. 

Intermolecular interactions included electrostatics, 
modeled using the effective charges approach of Gaba- 
doulline and Wade [44], and combined van der Waals 
and hydrophobic interactions described with the 
Lennard-Jones potential [42], with the well depth treated 
as an adjustable parameter. Hydrodynamic interactions 



were neglected. The authors investigated how different 
energy models affect the translational diffusion coeffi- 
cient of GFP. When only steric interactions (modeled 
with the repulsive part of the 6/12 Lennard-Jones poten- 
tial) were considered, the BD simulated translational dif- 
fusion coefficient of GFP was 3 to 6 times higher than 
the coefficient estimated experimentally. Even including 
electrostatics did not significantly reduce the discre- 
pancy between the experimental and computed values. 
Nevertheless, the authors obtained a close agreement 
with experiments after adding the short-range Lennard- 
Jones attraction and tuning its well-depth. The observa- 
tion that purely steric interactions are not able to 
account for the reduction in proteins' mobilities in vivo 
concur with previous computational studies [71]. 

While adjusting the short-range attractive interactions 
led to reproducing the in vivo GFP diffusion coefficient, 
it has been later pointed out [37] that the magnitude of 
attractive interactions can be tuned to set the diffusion 
coefficient to any chosen value. McGuffee and Elcock 
investigated also the character of translational diffusion 
[4] in cytoplasm with different interaction energy mod- 
els. When only the steric interactions were considered, 
for the three chosen proteins they observed normal dif- 
fusion at short times (ps), subdiffusion at moderate time 
scales (ns) and normal diffusion at long times (ps). 
When electrostatic and short-range attractions were 
included, sub-diffusive behavior was noticed even for 
short observation times. 

The approach to model diffusion in the E. coli cyto- 
plasm taken by Ando and Skolnick [37] is different than 
the one described above. These authors investigated the 
effect of macromolecular shapes on the diffusion in 
crowded, heterogeneous environment, the role of HI, and 
the interplay between hydrodynamic and non-specific 
attractive interactions. Their cytoplasm model consisted 
of 15 different types of molecules, each described with a 
set of spherical subunits positioned at the C a protein 
atoms and the P, C4', Nl, and N9 atoms of nucleic acids. 
The authors dealt with the anisotropy of macromolecular 
shapes using the rigid body BD algorithm [26] that 
describes molecules with precomputed, 6x6 diffusion 
tensors. Notably, this is the first BD study of a multimo- 
lecular system that takes into account the anisotropy of 
diffusion. To investigate the role of HI they also simu- 
lated the system of equivalent spheres corresponding to 
the cytoplasm model (all molecules in the system were 
represented as spherical particles and assigned hydrody- 
namic radii based on translational diffusion coefficients 
computed from their atomic structures [27]) with the far- 
field many-body HI interactions and lubrication forces 
treated with the method of Banchio and Brady [62] . Soft- 
sphere harmonic potential was used to model repulsive 
interactions between particles. Attractive interactions 
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were described by the Lennard-Jones potential function 
with a correction factor accounting for the roughness of 
molecular surfaces. Electrostatic interactions were mod- 
eled using the DLVO model. 

Based on simulations of molecular-shaped and equiva- 
lent sphere systems, the authors showed that the effects 
of macromolecular shapes on long-time translational 
diffusion are rather small in the studied concentration 

range of 250 - 350 ^S- . This is an important result that 
ml 

justifies the use of isotropic models to study transla- 
tional diffusion in the crowded systems. The influence 
of macromolecular shapes on rotational dynamics, 
which is an interesting issue as well, was not investi- 
gated. Their BD simulations with only steric repulsion 
between particles also resulted in overestimated transla- 
tional diffusion coefficients. However, when HI were 
switched on, the computed long-time translational diffu- 
sion coefficient of GFP was in good agreement with 
experimental value; this is another important result 
because this model did not include any tunable para- 
meters. The authors concluded that the macromolecular 
motions in vivo are likely to be dominated by steric 
interactions and hydrodynamics. Additionally, they 
observed that when non-specific attractive interactions 
between particles (either molecular-shaped or spherical) 
are included, the reduction of translational diffusion 
coefficients depends much more on molecular sizes 
in comparison with the simulations where HI were 
enabled. 

Both studies described above [37,43] are state of the art 
applications of BD in the field of macromolecular diffu- 
sion in biological systems, as they both take into account 
the heterogeneity of the environment and intermolecular 
interactions beyond simple excluded volume models. 
Both are also closely connected to experimental data. 

Concluding remarks 

This short review focuses on the current status of the 
particle-based Brownian dynamics technique and its 
applicability to model diffusive phenomena in biological 
environments. We described a few BD studies of three- 
dimensional macromolecular diffusion in aqueous cellu- 
lar compartments under molecular crowding conditions. 
However, there are various other BD applications, not 
presented here, that consider diffusive phenomena in 
cellular architectures, such as the cytoskeleton [72] or 
membranes [73]. 

We presented different modeling approaches, from 
simple excluded-volume models to most sophisticated 
atomic level descriptions of the in vivo-like systems. The 
complexity of intra and extracellular environments, their 
dynamics, and the number of different interactions that 
influence the diffusive behavior of macromolecules, 



make computational studies of crowding extremely diffi- 
cult. While it is desirable to construct models that treat 
biological systems at atomic resolutions and include 
sophisticated interactions [42,43], the recent model of 
the E. coli cytoplasm [37] shows that less elaborate, 
coarse-grained approaches should not be disregarded, at 
least as far as diffusive transport is concerned. 

Obviously, current BD algorithms need development. 
These include accurate and computationally efficient mod- 
els of direct, both specific and non-specific, intermolecular 
interactions as well as HI models that up to date have 
been largely neglected in the simulations of biological 
systems. 

BD simulations should be (at least ideally) able to repro- 
duce and predict the in vivo dynamics. However, the com- 
putational cost of sophisticated BD simulations may be 
too high, even with the growth of computing power and 
available modern technologies. Thus, the long-term goal 
should be to systematically develop efficient algorithms 
interfacing different modeling approaches - from atomistic 
to coarse-grained models of macromolecules, mesoscopic 
models of biological environments, and appropriate 
boundary conditions for different cell compartments. New 
approaches and models should be verified experimentally 
to establish their quality and predictive power. 
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